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1. Introduction 

The solution of the planar Dirichlet problem for the Helmholtz equation has served as 
one of the principal numerical setups for verifying and demonstrating the main ideas of 
quantum chaos [TJ [3], [4J, [5J, [6J, [7] . The reason is that this problem, usually referred to as 
the quantum billiard, is one of the simplest quantum Hamiltonian stationary problems 
for which a good quality spectra of highly excited eigenstates can be obtained for various, 
even non-integrable geometries. 

In literature there exists a good account of numerical techniques to tackle the 
problem [8j. Apart from many ingenious ideas developed for specific geometric setups, 
there are two competitive general purpose approaches: (i) Boundary integral method 
(reviewed in [9]) or (ii) Heller's plane wave decomposition method |2J. While the 
boundary integral method is really a general rigorously based method, the plane wave 
decomposition method is a rather heuristic approach which can fail in certain important 
general cases, e.g. of non-convex geometries [TD]- Still, there exists an extremely efficient 
implementation of plane wave decomposition method due to Vergini and Saraceno 
[12] , which makes this method an attractive option in spite of potential mathematical 
problems. The crucial new idea of Vergini and Saraceno was to look for a minimum 
of a boundary norm (an integrated norm of the wavefunction or the field amplitude 
along the boundary of a 2-dimensional domain) in terms of a solution of a generalized 
eigenvalue problem, the gain being that in a single computational step a number of 
accurate eigenvalues are obtained in a constant proportion to the dimension of the 
matrices scaling as O(k) where k = 27r/X is a referential wave-number, in contrast to 
the boundary integral method, where a number of matrix computations have to be 
performed in order to locate a single eigenvalue. In both cases each matrix computation 
is of the order 0(k 3 ). 

In this paper we propose to use a similar (scaling) idea in conjunction with 
the rigorous boundary integral method. We develop a completely general numerical 
boundary integral technique which produces a constant fraction of k eigenvalues per 
single matrix operation involving 0(k 3 ) scalar operations, thus being asymptotically 
orders of magnitude faster than traditional implementations of boundary integral 
method. 

The details of the method, which comprises the first part of the paper, are 
elaborated in Section 2. 

The second main idea of the paper is to apply our method in a chaotic quantum 
billiard whose phase space structure is not time-reversal invariant, in a sense that 
its chaotic phase space components are not mapped onto themselves upon the time- 
reversal operation. Such a situation can appear in generic convex billiards with generic 
KAM structure. Namely we study the billiard in the domain of non-convex and non- 
simply connected shape, the so-called Monza billiard. Our choice of the model system 
also represents a good benchmark for numerical methods since it is likely one of the 
most difficult types of the billiard shape that one can think of. The corresponding 
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classical billiard — being a curved closed corridor with parallel walls — possesses the 
property of unidirectional motion, namely the classical phase space separates into two 
disjoint ergodic and chaotic components of clockwise and counter-clockwise motions. 
The absence of any geometric (point) symmetries of the model and the fact that 
the time-reversal symmetry does not preserve the invariant ergodic components, has 
an interesting consequence, namely the existence of algebraically (in effective Planck 
constant ft) split quasi-degenerate pairs of energy levels. The effective Planck's constant 
of the system is defined as the ratio of the physical Planck's constant to a typical action 
of a system. In billiard systems, the effective Planck's constant % oc 1/k where k is the 
wave-number. The existence of quasi-degenerate pairs has a dramatic effect on energy 
level statistics, in particular since the level splitting is of the same order as the mean 
level spacing. The effect can be interpreted as a chaotic analogue of the Shnirelman 
peak [131 E] known for a long time for systems with non-time reversible KAM islands. 

Another and perhaps even more dramatic effect concerns the long range correlations 
among levels which behave according to the Gaussian Unitary Ensemble (GUE) of 
random matrices in spite of the fact that the system as a whole possesses time reversal 
symmetry. This is explained by the fact that even though the eigenfunctions of 
the Hamiltonian are strictly real, the doublets of eigenfunctions of the time-reversal 
operation T, namely T?^l,r = ^l,r> which become preserved in the classical limit even 
though they are not exactly eigenfunctions of the Hamiltonian, are complex and hence 
GUE should be used for modeling spectral statistics at large energy rangesJ3 This is a 
new effect, which to our knowledge, has not been predicted or known before. 

Detailed discussion of the Monza billiard, its spectral statistics and dynamics, is 
given in section 3. In section 4 we discuss our main result and conclude. 

2. Expanded boundary integral method 

The scaling method as proposed by Vergini and Saraceno [12] allows one to compute a 
number of states of a quantum billiard, or any other system described by the Helmholtz 
equation with the Dirichlet boundary condition, lying close to a chosen reference value 
of the wave-number k without losing any state in a chosen interval. One of its main 
advantages is that the matrices to be diagonalized are not of the dimension of the order 
of would be the case when using the usual diagonalization techniques, but only 

of the order oc k. While this fact makes it one of the most efficient approaches to solve 
for eigenvalues and eigenstates of quantum billiard problems, especially for high values 
of k, experience shows that its domain of applicability is rather limited. As shown by 
Gutkin in [10] , the plane wave decomposition method [2] on which the scaling method is 
typically based can in general only be applied to a convex billiard. The scaling method 
does not necessarily involve a plane wave decomposition and a different basis set can 
be used, which enables the method to solve problems inaccessible by the plane wave 
decomposition. One such example is the use of Bessel functions of fractional order to 

| For the general discussion of the role of antiunitary symmetries see [15] . 
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solve the mushroom shaped billiard problem (TTJ. The choice of the basis in such cases 
must, however, be specifically tailored to the system. 

The boundary integral method, see [9] for a review of the method, has a similar 
demand on matrix size as the scaling method and can in principle be applied to arbitrary 
shapes as it is based on the Green's theorem that shows how the wavefunction inside 
the billiard can be represented with its normal derivative on the boundary (in the case 
of Dirichlet boundary conditions). From this representation the integral equation for 
the normal derivative at the boundary u = -j^-ip follows as [9] 

u(s) = -2 / — G k (q(s), q(s')) u(s') ds>, (1) 
Jdn on 

where dQ denotes the boundary of the domain, s is the arc length parametrization of 
the boundary. The coordinates of the boundary are given as q(s) and -§^Gk (q, q') is 
the inward normal derivative of the Green's function for the unbound problem at wave- 
number k. Following reference [9] we may discretize equation |T| ona finite set of points 
as determined by their arc length parameters Sj, yielding 

[A(fc)u] i = X;^(*)«* = 0. (2) 

ikl 




+ "-^H? (A^)cos. 



(3) 



where Ui is the value of the normal derivative at the point q(si), Tij = \q(si) — q(sj)\ is 
the distance between points q(si) and q(sj), coscpij = n(s«) • (q(si) — q(sj))/rij gives the 
angle between the inward boundary normal n at the point q(si) and the distance vector 
between the two points and /tj is the boundary curvature at the point q(si), where 
positive curvature is taken for concave parts of the boundary. The arc-length of the 
boundary section centered at q(si) is denoted by The discretization is characterized 
by the parameter b = 2irN p /(kLiy), where N p is the number of discretization points and 
Lb is the total length of the billiard boundary. The parameter b measures the number 
of discretization points per wavelength. 

Equation [2] only has solutions for those values of k where the matrix A becomes 
singular. While in principle any Green's function could be used, the choice of the 
complex Hankel function of the first kind is necessary in order to avoid spurious solutions 
that occur because equation (TjQ) only then becomes a sufficient condition for determining 
the normal derivative. The usual approach towards solving this equation is based on 
solving for zeros of the Fredholm determinant. It therefore yields no more than a single 
level (if at all) per determinant calculation and is typically plagued with loss of levels and 
accuracy when close to a degeneracy. The latter problem has been resolved by Backer 
[9] by calculating not the determinant directly but the singular value decomposition 
(SVD) of the matrix A and following the behaviour of individual singular values. 

At this point we introduce a similar procedure to the one presented for the scaling 
method, namely trying to determine the solution close to a chosen reference value of 
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the wave-number fc by varying k = k + We may write the equation as a Taylor 
expansion, 



(6k) 



2 



A(fc ) + 5fcA'(fc ) + ^-^A"(A;o) + 



it 



0, (4) 



where A'(ko) and A"(ko) are the first and the second derivative of the matrix A with 
respect to k at the point ko, obtained by taking the derivatives of each matrix element. 
Let us expand 5k = e5ko + e 2 5k\ + . . . and u = Uq + eu\ + e 2 U2 + . . . in terms of 
the order of the formal perturbation parameter e, whose powers give the order of the 
terms with respect to the small parameter 5k and should be set to 1 in the final result. 
This also shows that the procedure can only be applied to those eigenvectors u whose 
wave-numbers lie close to the reference value k . By inserting these expressions into the 
equation (Tj0) and ordering the terms with respect to e we obtain (dropping the implied 
argument ko) 

Au + e5k A'u + tAui + 

+e 2 5k A'u 1 + e 2 5k 1 A'u + t 2 ^^A"u + e 2 Au 2 + O (e 3 ) = 0. (5) 

We will now show that taking the first two terms of the above equation, 

A(k )u = (-eSko) A'(k )u , (6) 

gives a consistent starting point for solving the perturbative problem. We may do so 
as there is a degree of freedom involved as to how to define the terms Uo and Ui in 
the perturbation series, and the equation (jSJ) therefore defines u . [jj This equation also 
shows that, while A and Uq are themselves of the zero order in the e expansion, their 
product is of the order 0(e). This is due to the fact that we are only slightly shifted in 
terms of wave-number k from the exact solution of the equation (j2J), where this product 
is exactly 0. 

Looking again at all the terms of the equation (jSJ) with up to the linear order in e 
and using the equation (jHJ), we obtain 

eA(k ) Ul + O (e 2 ) = 0. (7) 

As the matrix A is not singular unless we already chose ko to be the solution of 
our problem, the above can be satisfied only for Ui = 0. We can therefore write 
u = u + O{e 2 ). 

We may improve the accuracy of the levels by multiplying the equation (jSJ) with 
the adjoint of the left eigenvector vq, defined as 

[A(k )fv = (-e5k )[A'(k )fv . (8) 

Due to the equations (El) and (jHJ), the first four terms are exactly eliminated in pairs. 
Although A and vo are both of the zero order in terms of e, the product VqA = 0(e) 

§ In analogy to the standard problem of (quasi-)degenerate perturbation theory of quantum mechanics, 
we do not yet know even the starting eigenvectors Uq of the problem, and taking only the leading term 
is not sufficient as they can only be determined by the splitting due to the perturbation. 
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is of a higher order, as can be seen from equation (jSJ). This means that the last term 
in the equation (jSJ) when multiplied by the left eigenvector also becomes of the order 
0(e 3 ). Only two terms remain, giving 

, h (Sk ) 2 v ^A"(k )u 
1_ 2 v^A>(k )u [J) 
and therefore k = k + e5k + e 2 5k 1 + O (e 3 ). 

By calculating the equations (|SJ) and © for various values of k , the chosen 
spacing between the values depending on the desired accuracy, we may obtain all 
levels of a system within some chosen interval. As the average density of states varies 
approximately as oc k , by taking the above error estimate into account, the spacing 
between various reference wave-numbers ko at which the individual calculations are 
performed must vary as Ak oc ko {1/3) if the error of calculating an individual level is to 
remain a constant percentage of the mean level spacing. The method itself is applicable 
even for calculating the ground state of the system. 

For simply connected domains the above procedure is found to be robust with 
respect to the shape of the boundary, even in the case of arbitrary corners. The domains 
with holes require an additional step as spurious solutions are found in the spectrum. 
These solutions appear because the transpose of the matrix in equation (j2J), which is 
here used to solve the interior Dirichlet problem, gives the solution to the exterior 
Neumann problem as well. As the holes of the domain represent the exterior of our 
problem and are themselves compact, the left eigenvectors of our spurious solutions 
are simply the boundary values of the wavefunctions of the Neumann problem within 
the holes. We may therefore eliminate the spurious solutions by solving the Neumann 
problem for each of the holes by performing the same computation for each hole, where 
the boundary taken is now just the boundary of the hole in question, and then discard 
the levels obtained for each hole from the spectrum of the total problem. 



3. Monza billiard 



The method was applied to a billiard that was named the Monza billiard due to its 
similarity with the famous Italian racetrack. It comes from a family of unidirectional 
billiards, as defined in [16J. These billiards are shaped as channels with parallel walls 
and have the property that the classical trajectories going in one direction along the 
channel cannot reverse this direction. The phase space of such systems is therefore 
split into two disjoint regions, and the motion within each separate region can be fully 
chaotic. The two regions are separated by a one dimensional family of marginally stable 
bouncing ball orbits. The shape of the system as shown in figure [1] is one of the simplest 
nontrivial closed shapes without any geometric symmetry that belongs to this class of 
systems. For all the subsequent examples shown, the system parameters as defined in 
the figure were chosen to be q = 1/12, a = 1/2, b = 1/3, r = 1/3, a = 1 . The billiard 
was numerically tested to be classically fully chaotic and ergodic within each of the two 
invariant phase space components [T7] . 
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Figure 1. The definition of the Monza family of billiards. The external half-circles 
have the radius of 1, whereas for the internal ones the radius is q. The lengths of the 
straight segments are a and b. The straight segments are joined by circular arcs with 
the angle of a, and the shortest of those arcs has the radius r. The variables w 6 [0, 1] 
and z E [0, 1) (cyclic) as shown are used to parametrize the billiard. The above shape 
corresponds to the parameters chosen in our computation. 



We used the expanded boundary integral method to compute all levels of the Monza 
billiard up to the wave-number k = 70. We used the Weyl formula [18J to test that we 
did not lose any levels in the chosen interval. The discretization parameter was chosen 
to be b = 12. The spacing of the reference wave-number ko between individual runs of 

— 1/3 

the method was chosen to be Ako = 0.05A; 

The principle of uniform semiclassical condensation (PUSC, [19], [Tj, [20] ) states that 
in the semiclassical limit of effective h — > the individual eigenstates of a system 
correspond to invariant objects in the classical phase space. One would perhaps naively 
expect the eigenstates of a Monza billiard to correspond to either one of the two chaotic 
components. This, however, can not be the case since states corresponding to either 
chaotic component would necessarily have a nonzero probability current [15] . The 
system, on the other hand, possesses the time reversal symmetry and therefore its 
non-degenerate eigenfunctions must be fully real and as such can have no probability 
current. 

This seeming contradiction is resolved by noting that in the Monza billiard 
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most eigenenergies are to be found in nearly degenerate doublets. Only the states 
corresponding to the bouncing ball modes are singlets, with their relative measure 
vanishing in the semi classical limit, as they correspond to a classical invariant separatrix 
between the two chaotic components which itself is the family of bouncing ball orbits 
corresponding to the particle bouncing perpendicularly between the two walls and is 
of measure in the classical phase space. The doublets are a chaotic manifestation of 
the Shnirelman peak in the level spacing distribution [13J that is generally found in all 
systems with a time reversal symmetry but having no point symmetries. Typically 
in such systems, the Shnirelman peak is due to the states that correspond to the 
classical invariant tori which do not map back onto themselves upon the operation 
of time reversal. There is usually exponentially small tunneling (in terms of effective h) 
between the tori that are connected via the time reversal operation and this is reflected 
in an exponentially small energy splitting between the states that correspond to the 
torus pair. In the Monza billiard, however, the tunneling does not occur between tori 
but between two chaotic components that are separated by a bouncing ball manifold of 
measure in phase space. This means that the expected tunneling effects will not be 
exponentially small but will rather scale as a power law in terms of effective Ti as the two 
chaotic components are a distance apart. In fact, it appears that the average splitting 
of the pair remains constant as a fraction of the mean level spacing. Heuristically, this 
can be explained by noting that the tunneling amplitude between the states localized in 
classically invariant chaotic phase space components can be estimated in terms of the 
phase space overlap of the corresponding Wigner functions, 



which can be semiclassically estimated to scale as a linear function of an effective Planck 
constant. 

To demonstrate this, let us fix an arbitrary point in the billiard q and locally 
represent the wavefunction as a combination of random plane waves, 



where p v = p (cos(ip+j3(q)),sm(tp+j3(q))). Due to the unidirectionality, only the waves 
in one half of the wave- vector space are taken, where /3(q) denotes the polar angle of 
the bouncing ball trajectory going through the position q. We assume a distribution of 
the stochastic variable / such that 




(10) 






with p = —t where A is the area of the billiard and is the Heaviside step function. 
Inserting the expression ffTT]) into the definition of the Wigner function 




W(q,p) = 2 / ^ x - X / 2 )HQ + x/2) exp(-ip ■ x/h) (13) 
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•/»iV) exp (i {{p v + zv)/2 - p) ■ x/%) . (14) 
We define the 2D domain 

V(q) = {x; (q + x/2 E ft) A (q - x/2 E fl)}. (15) 
The expected value of such a Wigner function is then 

(W L (q, p)) = -—72 / d 2 x f d V p exp (i (p - p) • aj/n) . (16) 

(27T^,) JO 

Performing the integration with respect to x the exponential turns into a wide delta 
function 

5 a (p) = — ^ [ d 2 ' 

(27lh) 2 Mq) 

such that 



MP) = 77TTT2 / ^ a;exp(ip-a;/fc), (17) 
(27m) -^(g) 



/*7T 

(Wl(9, p)) = / dip p 5 a (p v - p) 
Jo 



The typical width o of the delta function is nonzero due to the finite size of the system. 
Its value is o oc h/£ where i is some linear dimension of the system independent of h. 
To get an estimate of the overlap in the equation (TTOT) we also need to introduce the 
Wigner function of the state traveling in the opposite direction 

(W R (q,p)) = !^ d V p 5 a ( Pip -P) (19) 

J TT 

with the integration boundaries complementing those in the expression (fT8]) . If we 
neglect the correlations we may use the expected values for the Wigner functions in the 
equation ([TO]) . The expected values of the two Wigner functions then overlap due to the 
width a, and the overlap itself is proportional to that width, thus demonstrating the 
statement (JTUJ). 

The spectrum can therefore be thought of as a composition of two nearly identical 
sequences, where their difference remains significant on the scale of the mean level 
spacing. It is important to stress that this near degeneracy can not be removed by any 
means such as desymmetrizing of the billiard as the system chosen does not possess any 
point symmetries. From this we can deduce that the unidirectional family of billiards 
will generally exhibit a non-universal spectral statistics. In the level spacing distribution 
P(S), which is the distribution of spacings S between consecutive levels in an unfolded 
spectrum where the mean level spacing is equal to one, we expect two main contributions, 
namely one from the spacings within individual pairs and another from the spacings 
between consecutive pairs. The first contribution is a widened delta distribution with 
the weight 1/2 close to the spacings 5 = 0, and the other contribution would be expected 
to be a stretched GOE (Gaussian Orthogonal Ensemble [2TJ) contribution. The widening 
of the second contribution by a factor of two is due to the unfolding procedure which 
requires that the mean level spacing is equal to one. In a sequence of pairs, the average 
spacing between centers of different pairs is therefore twice the mean level spacing. 

|| Note that this estimate is independent of geometrical dimension of the system, i.e. it should be the 
same for 3D "tube" billiards. 
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The cumulative level spacing distribution 

W(S) = [ P(S')dS', (20) 



for the system with the parameters q = 1/12, a = 1/2, b = 1/3, r = 1/3, a = 1 is given 
in figure [2j The number of levels used in the figure was iV = 2403. For comparison, 
the same statistics using only the last 500 levels is shown (dashed red) as well in order 
to demonstrate that the level spacing distribution indeed does not change with energy 
(or h). Apart from the initial contribution that is mostly due to the spacings within 
individual pairs, the level spacing distribution then follows the stretched and GOE 
prediction 

rS/2 

1 + / P G oE(S')dS' , 
Jo 



W G oe(S) = \ 



(21) 



where 

7T „ / 7T 



Pgoe(S) = exp (-^S 2 ) (22) 



is the Wigner surmise [21]. This prediction holds well for small and intermediate 
spacings, S < 2 (shown dotted green), yet at larger spacings the distribution starts to 
deviate from the GOE prediction but rather approaches the (stretched) GUE prediction 
(dot-dashed blue). This can be obtained from equation ( 12TI) but replacing Pgoe with 

P GU E(S) = - 2 S 2 exp(--S 2 ). (23) 

While not exactly following this prediction, a clear trend does emerge for the W(S) to 
move from the GOE to the GUE behaviour at larger spacings. 

The longer range spectral statistics are investigated by using the spectral rigidity 
[211. which is defined as 



A 3 (L) = (min - ^ [N(x) - (Ax + B)f dxj , (24) 

where N(E) is the spectral staircase function whose value increases by one at each 
(unfolded) energy level. The rigidity measures the average deviations of the spectral 
staircase function from the best fitting straight line of length L. It is given for the 
spectrum of the Monza billiard in the figure [3J We are comparing the numerical results 
to the A GO (u)e predictions for the random matrix ensembles as defined in [21]. Since 
most levels come in pairs, the predictions need to be scaled as 

&go(u)e{L) = AA go{U )e(L/2). (25) 

The stretching in L is necessary because the average level spacing between level pairs 
is twice the average level spacing. The factor of 4 occurs since each step of the spectral 
staircase function is twice as large for each pair as it would be for an individual level, 
and the A 3 statistics is quadratic in spectral staircase fluctuations. 

The results for the A 3 statistics of the Monza billiard are given in figure [3j At low 
values of L, the results of both predictions as well as the numerical results agree. The 
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S 



Figure 2. The integrated level spacing distribution W(S) for the Monza billiard. 
The numerical results for the 2403 lowest levels are shown as the full black line. The 
results for the highest 500 levels of the same sequence arc shown as a red dashed line. 
The stretched and shifted GOE (dotted green line) and GUE (dot-dashed blue line) 
predictions are shown as well. In the inset the difference of the GUE and GOE W(S) 
distributions to the numerical results is shown. 



numerics undershoots both predictions because the pairs are not exactly degenerate, 
with the small spacings within the pairs reducing the fluctuations of the spectral 
staircase function in comparison to the assumed exact degeneracy in the GOE and 
GUE predictions. At larger distances L, the numerics clearly follows the GUE rather 
than the GOE prediction. At the distances larger than about L = 10 the numerical 
A 3 statistics starts to increase faster than any of the predictions. This effect could be 
expected because of the inherent presence of the large bouncing ball mode contribution 
to the spectrum [22]. Excluding this effect, however, the behaviour of the spectrum at 
larger distances appears to follow the GUE prediction. 

To understand such behaviour we need to look into the properties of state pairs. 
In figure 0] we show two states corresponding to a neighbouring pair of states ip a and 
if}},. Both states are of course fully real, with the phase switching between (shown red) 
and 7r (shown cyan), but a linear combination ipL,R = {if> a ± iipb)/y/2 of two real states 
can be shown to have the largest current amplitude. In our case, the direction of the 
phase change, which itself corresponds to the direction of the current probability, shows 
that this current flow is clearly seen to be predominantly in a single direction. The 
complex conjugate of the function is orthogonal to the function itself and corresponds 
to a current in the opposite direction. The current behaviour is even more clearly 
exposed by looking at the Husimi plots of the ID wavefunction cross section along the 
coordinate z taken at w = 1/2 (see figured] for coordinate definitions). While for the 
two eigenstates the Husimi plots are nearly identical, examining the Husimi plot for 
the complex combination of the two states clearly yields a state that corresponds to a 
current in a single direction. This is typical for most pairs of states and, while PUSC 
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Figure 3. The As(L) spectral rigidity for the lowest 2403 levels of Monza billiard is 
shown as a full black line. The scaled GOE (dotted green) and GUE (dashed blue) 
predictions are given as well. 



can not be strictly applied to individual eigenstates, taking a nearly degenerate pair and 
performing the above complex rotation in this subspace yields states that support the 
PUSC conjecture^ 

The nature of the spectrum naturally poses the question of the dynamics of such 
a quantum system. Classically, a particle traversing the billiard in one direction will 
keep that direction permanently Starting with a quantum state such as shown in the 
bottom part of the figure HI as this state is not an eigenstate of the system, it will 
exactly reverse its character every t = nh/AE where AE is the energy splitting of the 
pair of the two eigenstates forming the initial state. Let us now define the operator for 
the current along the billiard direction as 

J = -2ie z ■ V, (26) 

where e z is the unit vector in the direction of the z coordinate in the billiard (see figure 
[T]). The initial state is chosen to be a Gaussian wave-packet 

V>(r, t = 0) = — exp (--^j- + j (27) 

with p Q = (10,0), r = (0,1/2) and a = 1/5 with the coordinate frame as shown in 
figure HI We expand this initial state in the basis of all eigenstates up to the wave- vector 
k = 20, where these states were obtained using the boundary discretization parameter 
b = 18. The evolution equation is taken to be the dimensionless Schrodinger equation 

id t ip = -Aip. (28) 

The expected value of the current operator J as a function of time is shown in figure 
|5j As can be seen, for a very short time the initial wave-packet has a large current 

A similar analysis applies to (almost) degenerate tori in a KAM-like scenario. 



Expanded boundary integral method 



13 




-2 -1.5 -1 -0.5 0.5 1 15 o 0.Z 0.4 0j6 0.8 1 




-2 -1.5 -1 -0.5 0.5 I \5 $ 0.1 04 0j6 OS 1 



Figure 4. States corresponding to the nearby pair of states with k a = 60.0725 (top) 
and k b — 60.0740 (middle), along with their complex linear superposition (bottom). 
The Husimi distributions for the wavefunction cross-section along the middle line of 
the Monza billiard (w = 1/2, z G [0, 1), see figure [T]) are shown to the right of each 
corresponding wavefunction. The darkness of the wavefunction corresponds to the 
probability density, whereas the hue going circularly from red to green to blue to red 
again represents the change of the phase of the wavefunction. In the Husimi plots 
there are 10 contours of constant value shown, starting from to the maximum value, 
where on the x axis the coordinate that runs along the length of the Monza billiard z 
is given, while on the y axis the momentum along this direction is shown in the units 
of the maximum attainable momentum at the corresponding classical energy. 

that up to until the time t = 1 relaxes to the value vq2/ti (shown dashed green), where 
v o = 20 is the classically expected velocity for the wave-packet with the chosen average 
energy. Classically, the current is expected to remain at this plateau, but the quantum 
system experiences a decrease of current that actually turns to fluctuating reversals of 
the current as it does not stabilize at 0. This behaviour remains even for times much 
larger than the ones shown in the figure. The fluctuations are seen to be a fairly large 
proportion of the initially chosen current. 

The stills from the video showing the behaviour of the wavefunction for the times 
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Figure 5. The expected value of the current operator as a function of time is shown as 
the full black line. The inset shows the same plot magnified for short times. The dotted 
blue line shows the value while the dashed green line gives the classically expected 
current for the whole chaotic component. The vertical dotted red line indicates the 
Heisenberg time tu — where AE is the mean level spacing. 

up to t = 1 are shown in the figure EJ A slowed down representation for the times 
up to t = 0.3 is shown in the figure CD These show the evolution of the wave-packet 
as it spreads out up to the point where it almost uniformly covers the billiard and 
therefore corresponds to the wave-packet as spread over the whole unidirectional chaotic 
component. The phase of the wavefunction clearly indicates that the direction of the 
quantum current is indeed directed in a single direction along the billiard for all the 
times shown. In figure [8] we show a still from the video representing the evolution of 
the wavefunction at the time t = 175 until the time t = 175.3, where the current as 
shown in the figure [5] experiences its strongest reversal. From the phase representation 
the direction of the current is not easily read. The figure [5], however, shows that the 
global unidirectional current contribution is quite sizable but the local fluctuations seem 
to obscure it. 

4. Conclusions 

The goal of the present paper was twofold: 

(i) Presentation of a novel efficient method for numerical solution of the 2D 
Helmholtz equation with Dirichlet boundary condition (e.g. a quantum billiard) which 
is obtained by combining rigorously founded boundary integral method with the ideas 
of Vergini and Saraceno [12J of expanding the boundary norm as a function of energy. 

(ii) Application of the new technique to compute the energy spectra, their statistical 
properties, and time evolution in a quantum billiard with a chaotic classical limit whose 
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Figure 6. The snapshots of the video ( |GIF animation, 8.0MB[ ) showing the wave- 
packet evolution (defined in text) for the times t = (top left), 1/4 (top right), 1/2 
(bottom left), 3/4 (bottom right). The representation is the same as in figure |H 




Figure 7. The snapshots of the video ( |GIF animation, 6.2MB[ ) showing the wave- 
packet evolution (defined in text) for the times t — (top left), 0.09 (top right), 0.18 
(bottom left), 0.27 (bottom right). The representation is the same as in figure H) 
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Figure 8. The snapshot of the video ( |GIF animation, 8.9MB[ ) showing the wave- 
packet evolution (defined in text) for the time t = 175. The representation is the same 
as in figured] 

phase space structure is not time reversal invariant. Here we mean that the system 
itself is (globally) time-reversal invariant but different disjoint chaotic components are 
not mapped onto themselves upon time reversal operation. Such situation is possible 
in the so-called Monza billiard due to unidirectionality of the classical motion [16] . We 
show the presence of energy level statistics converging to a well defined distribution 
in the semiclassical limit, which is non- universal and exhibits a Shirelman-type peak 
at small energy ranges and GUE fluctuations at large energy ranges. We have also 
computed the time evolution of the expectation value of the tangential momentum 
(or particle current) which exhibits interesting irregular oscillations due to algebraic 
quantum tunneling phenomena. The time-scale of the current reversal and consequent 
oscillations is a constant factor (typically much larger than unity) of the Heisenberg 
time. The results have been qualitatively explained in terms of a heuristic semiclassical 
picture, namely the principle of uniform semiclassical condensation combined with the 
nontrivial effects due to time reversal symmetry. 

The phenomenon of GUE-like fluctuations in a time-reversal invariant system, 
which we discovered at large energy ranges, is similar to the behaviour observed by 
Leyvraz, Schmit and Seligman [23] in systems with point symmetries lacking real 
representations. Indeed the two phenomena can be understood to have some formal 
similarities. In future explorations it would be interesting also to consider systems 
which live in both situations simultaneously, having classically unidirectional motion 
and point symmetry with complex representations. 

As the quantum billiard models can nowadays be easily realized in various 
experimental setups (e.g. quantum dots, microwave cavities, optics etc) we expect that 
predicted effects should be experimentally observable. 
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Appendix A. Merging level sequences 



To obtain the whole spectrum from the individual runs at various reference values k Q 
a method must be chosen to merge all the individual level sequences into a single list. 
While in most applications a simple trimming of individual sequences so that they fit 
with their neighbours on the energy scale is sufficient, there always exists a possibility 
that, if a level is very close to the trimming boundary of a subsequence, it may be 
accidentally removed, or doubled if, due to the numerical inaccuracy, it happens to be 
present just barely within the neighbouring list as well. To merge two subsequences 
and the first calculated at a reference value ko and the second at ko + Ako, we 

first define a weight function p(x), where x(k) = Ak ~ 1 (k — ko)- The function is defined 
such that p(x < 0) = pix > 1) = 0, p(l/2) = 1 and that its derivative is bounded such 
that \p'ekAko\ <C 1, where e& is the typical numerical error of an individual level. We 
used the function p(x) = 4x(l — x). First we calculate 



c — - 



which gives the estimate of the (weighted) number of levels in the interval (k , ko + Ak ). 
If this weighed estimate is smaller than, say 1/2, we know there are no levels in the 
neighbourhood of the middle point of the interval and we may just merge the two 
subsequences by trimming them in the middle and joining them. Otherwise we trim the 
first sequence by discarding all the levels that are higher than the middle point of the 
interval and then vary the trimming point k c for the second sequence in such a way that 

a'(k c )= £ p(x(k i ))+ £ p(x(k' j )) (A.2) 

_{ki<k +^} { k 'j> k -} 

differs as little as possible from a. One should also always demand that k c > k 
(preferably even k c > k + AMsl to avoid small contributions to a' from the edges of 
the interval), otherwise possible levels outside the chosen interval may be inadvertently 
added. 
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